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Abstract 

We consider general multi-species models of reaction diffusion processes and 
obtain a set of constraints on the rates which give rise to closed systems of 
equations for correlation functions. Our results are valid in any dimension and 
on any type of lattice. We also show that under these conditions the evolution 
equations for two point functions at different times are also closed. As an example 
we introduce a class of two species models which may be useful for the description 
of voting processes or the spreading of epidemics. 
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1 Introduction 



An interesting class of non-equilibrium problems with a rich dynamical behaviour and 
a vast area for applications are stochastic reaction-diffusion systems (see 0, ||, H, 
j|] and references therein). These are the processes in which one or several species of 
particles hop randomly on a lattice, and interact in various possible ways with each 
other. In the one species case, it is convenient to denote a particle by the symbol 1 
and a vacant site (a hole) by the symbol 0. Then a simple hopping is represented 
by 1 + 0— ► + 1. In addition to exclusion which means that no two particles can 
occupy a single site, the possible interactions include pair creation (0 + 0— >1 + 1), 
pair annihilation (1 + 1— »0 + 0), coagulation ( 1 + 1 — > 1 + 0), de-coagulation, 
(0 + 1— > 1 + 1), birth ( — > 1) and death (1 — > 0) processes. Obviously the variety 
of elementary processes rapidly increases with the number of species. 
In general such lattice systems are difficult to treat by rigorous analytical means and 
correspondingly, considering the vast amount of such models, relatively few exact re- 
sults are known. 

Over the past few years the application of operator formalism to these stochastic 
processes and their mapping to quantum spin systems and their generalizations has 
turned to be quite fruitful. In view of this correspondence, many of the techniques 
of quantum spin systems such as free fermion techniques, Bethe ansatz and related 
algebraic techniques || |(| |?J ^, ^, |ll| , have been used in the investigation of reac- 
tion diffusion systems, specially in one dimensional lattice systems. Also by using the 



operator formalism, some of the old techniques such as the matrix product ansatz [12 



have been put to very fruitful use, in solution of one dimensional stochastic systems 
[ ]D1 ] (see jl], ^, 14 1 and references therein). Almost all of the above methods have one 
limitation, they are restricted to one dimensional lattices, specially if we are interested 
in exact solutions. 

A common feature of any model of interacting particles and indeed the main source 
of difficulty in obtaining exact solutions is that the equations of motion of correla- 
tion functions form an infinite hierarchy, that is, the equation of motion of n— point 
functions include (n + 1)— point functions and in general higher correlation functions. 
One can truncate the hierarchy at a level by various kinds of approximations, the 
simplest and the most common method is to break the hierarchy at the first level by 
the mean field approximation. However in some models and in low space dimensions, 
the amount of diffusive mixing may not be enough to warrant such an approximation. 
Some of these models may also have long relaxation times so that their simulation 
may be difficult and time consuming. For these models exact solutions are highly 
desirable. 

Interestingly enough, there are models in which the hierarchy of equations of correla- 
tion functions automatically breaks at every level, hence the possibility of obtaining 
exact solutions. For these models the equations of motion of n— point functions de- 
pend only on k— point functions with k < n. We should stress that while a powerful 
technique like matrix product ansatz or its dynamical version [15] only transform the 
problem to another equally difficult and yet convenient problem, that is, finding the 
representations and calculating matrix elements of strings of operators of an alge- 
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bra [13, 16, 17], for these models the simplification appears to be genuinely effective. 
Moreover this property is independent of the dimension and the geometry of the lat- 
tice. In view of its generality this is a great simplification and deserves to be pursued 
further from various directions. 

The observation of this phenomenon in some models like symmetric exclusion and 

partial exclusion processes @, m, m in ied g - m - schutz m to raise the question 

of classification of such models, i.e., a general criterion on the reaction and diffusion 
rates such that the resulting equations for the correlation functions decouple. By 
considering the one species processes and the particle density correlation functions 
(n(x)), he found that from among the 12— parameter family of single species reaction 
diffusion systems, a 10— parameter family fall within this class. 
Since then this question has been pursued further. For example one can consider the 
hole-density correlation functions (1 — n(x)) (the so called empty interval method) 

or even more general functions like 



with two [23, 24] or three site interactions 

, |26f| . These correlation functions lead to different sets of constraints 



(a + bn{x)) p2| 
on the rates 



1.1 The aim and the results of the paper 

The aim of the present paper is to investigate this question for the general multi- 
species case. That is, we assume that there are p + 1 species of particles labeled 
as type 0, 1, 2, ■■■p which hop and interact in a lattice of arbitrary geometry. We 
interpret particles of type as holes and other particles as real particles. We assume 
two body interactions between the particles, and obtain the general condition on the 
rates, so that the correlation functions of densities of real particles, i.e., particles of 
type 1, 2, • • • p decouple from the correlation functions of higher levels in the hierarchy. 
We show that under this condition the equation of motion of two-time two-point cor- 
relation functions are also closed. These types of correlation functions are important 
in the analysis of voting processes [27]. 



In view of the arbitrary number of species in our model, one will have much more 
freedom to find an exactly solvable model for description of physical phenomena, spe- 
cially in the area of chemical kinetics, where one usually has more than one species 
of particles. 

Since in our analysis and hence in our results, there is an asymmetry between a parti- 
cle of type and other particles, it is important to keep in mind that the interpretation 
of particle as a hole is not essential. Therefore in adapting a model from the class 
discussed below to a situation of physical interest, one can change this interpretation 
and take any other particle of type 1,2, • • -p to stand for a hole and particle for a 
real particle. 

For general p— species models, the number of independent rates is (jp + l) 4 — (p + l) 2 . 
We show that the family of exactly solvable models (in the above sense) live on an 
N p dimensional manifold (hyperplane) where N p := (p+ l) 4 — (p+ l) 2 — 2p 3 . For the 
one species case we find N\ = 16 — 4 — 2 = 10, in accord with the results of Schutz 
plj|. For the two and three species cases we have respectively N2 = 56 and A^3 = 186, 
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where the number of independent rates for these models are originally 72 and 240 
respectively. For an elaboration on this see the remark 3 in the text. 
As an example and for concreteness we study a 2— species family and further con- 
strain it with certain extra symmetries. In this class we can find models suitable for 
the description of the spread of an epidemic, the exchange of ideas and votes, and the 
spreading of news or rumor. 

We set up the solution of one point functions for these models and investigate to a 
certain extent the properties of these solutions. 

The structure of this paper is as follows: In section 2 we introduce our notations and 
conventions. In section 3 we obtain the general conditions on the rates. In section 4 
we consider the case p = 2, and by imposing further symmetry requirement on these 
models we introduce a class of solvable two-species models. In section 5 we set up 
the general solution of the one point functions for this model. We conclude the paper 
with a discussion. 



2 Notations and conventions 

Throughout the paper we will adhere to notations which we will collect here for 
convenience. 

For ease of notations we will consider a one dimensional lattice which may be infinite 
or periodic. However all our results are valid also on lattices of arbitrary shape and 
arbitrary dimension. We will also work explicitly with nearest neighbor interactions, 
although again our results are valid for arbitrary range of interactions. These facts 
have already been shown in and are also easily verified by reviewing our method 
of proof, in the sense that no part of the reasoning depends on the underlying lattice 
or the range of interactions. 

We denote the points of the lattice by Latin letters from the end of the alphabet 
x — 1, x, x + 1, • • •. In a finite lattice we number the sites from 1 to L. To each site x 
of the lattice, we assign a random variable t(x) which can take p + 1 values 0, 1, • • • p. 
We denote the values of this random variable when it takes all the values including 
possibly the value 0, by Greek letters a, (3, \i, ■ ■ •, and when it takes only the values 
different from 0, by Latin letters, from the middle of the alphabet, like i,j, k, 1-. Thus 
(rij(x)) := {8 T r x \i) denotes the average density of particles of type i at site x, or the 
probability of site x being occupied by a particle of type i, and (no(x)) := (8 T ( x )fi) 
denotes the probability of this site being empty. It thus follows that 

£ < n,{x) >= 1. (1) 
fi=0 

We assume that two particles of type a and (3 on two adjacent sites may transform 
stochastically to particles of type /i and v with rate -R^' This is written as: 

a,P — > fi, v with rate R^g- (^) 
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Since a Greek index includes also the value 0, which we interpret as a vacant site, the 
above transformations include all the possible processes conceivable for all types of 
particles on the two sites. For example R®q is the hopping rate of a particle of type i 
and Rqq is the rate of creation of a particle of type k from the vacuum and Rfj is the 
rate with which two particles of type i and j interact or coagulate to form a particle 
of type k. 

We will use the operator formalism for Markov processes. This formalism is well 
known by now, since in the past few years it has been used extensively for the analysis 
of reaction diffusion processes, specially in one dimension. 

In this formalism we should assign a complex p + 1 dimensional Hilbert space C p+1 
to each site, with ortho-normal basis states 

|0>,|l),--», W) = &w (3) 

The Hilbert space of the whole lattice is the tensor product of all the local Hilbert 
spaces of the sites. At any given time t, each configuration of the lattice is given 
by the values of the random variables of all sites r(l),r(2), • • • r(L). Such a config- 
uration occurs with probability P(r(l),r(2), • • • t(L); t) = (r(l),r(2), • • • r(L)\P(t)). 
The state vector \P(t)) determines all the probabilities and its evolution is governed 
by a Hamiltonian derived from the rates: 

j t \P} = »\P) (4) 

For convenience we have absorbed the minus sign, which is usually included in this 
equation, into the definition of the Hamiltonian. The probabilities are normalized by 
requiring that 

£ (r(l),r(2),...r(L)|P) = l (5) 

r(l),r(2),-r(L) 

This equation can be rewritten as 

(S\P) ■- (s[® L |P) = 1 (6) 

where (s| is defined as: 

<'|:=f>l (7) 

Note that the Bra state (S\ is the sum of all the possible configurations of the system. 
From the master equation (|4|) and the property {S\H = (as required by conservation 
of probability), one obtains the Heisenb erg-like equation of motion for the average of 
any time independent observable O: 

j t (0(t)) = (S\[6,H]\P(t)) (8) 

where O is an appropriately chosen operator whose matrix element gives the average 
of the observable, i.e: 

<0(t)> = (S\6\P(t)). (9) 



4 



Note that the time dependence of the average comes from the evolution of probabili- 
ties, therefore a better notation for the average of a time independent observable will 
be (0)(t). However we use the notation (0(t)), since it will be convenient when we 
consider two point functions at un-equal times. Consider a completely general Markov 
system whose configurations are labeled by C and a time independent observable O 
of the configurations. We have: 



(0(t 2 )0(ti)> := ]T 0(C 2 )0(Ci)P(C 2 ,t 2 ;Ci,ti) 

C?2,Cl 



(10) 



In the operator formalism, it is easily shown that this two point function is represented 
as: 



(0(t 2 )0(tx)) := (S\0(t 2 )0(h)\P(p)) 



(11) 



where 0(t) is the Heisenberg-like operator, 0(t) := e Oe . The above cor- 
respondence is also true for different kinds of operators. In particular we have: 
(O(t)O(O)) := (S\e- tH 6e tH 6\P{0)}. We then have: 



d_ 
~dt 



(O(t)O(O)) = (S\[O(t),H]O\P(0)) = {S\[O,H]e tH O\P(0)) 



(12) 



where we have used (S\e~ tH = (S\. 

Returning to our model, the local operators E a p := \a){/3\ act on the states of 
a site as: E a ^\(j) = 5p^\a). Of particular interest are the diagonal operators 
Eqq, En, ■ • • E pp which act as number operators for the holes, particles of type 1, 
to particles of type p. Hereafter we abbreviate these diagonal operators and show 
them only by one index instead of two, i.e: E„ stands for 



En 



( 1 



V 



o 



o 7 



E n 



( o 



V 



o 



\ 



i / 



(13) 



The above operators have the commutation relations: 

[E a p,E^ u ] = 5fj^E au — 5 l/a Epp. (14) 
Moreover the following property of these operators is also important in the sequel: 

{s\E aP = {s\E pp =: { 8 \E P , (15) 
where we have abbreviated Epp as Ep. 

The Hamiltonian describing the processes @) is constructed as a sum of local Hamil- 
tonians acting on adjacent sites: 



(16) 
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where the operator h(x, x + 1) means that the operator h acts only nontrivially on 
sites x and x+1. The operator h is constructed from local operators as follows, where 
we will use hereafter the Einstein summation convention over Greek indices: 

h = R^(E m ®E v p-E a ®E l3 ) (17) 

The conservation of probability constrains the rates to satisfy the relation 

$>$ = V a, (18) 

Finally we need two matrices constructed from the matrix of rates which will prove 
useful later, and we prefer to introduce them here for convenience. Out of the matrix 
of rates, we form a set of matrices R l and 5* defined as follows: 

(R%? ■= £iC/3 (19) 

V 

(S%,p := Y, R «P (2°) 

V 

The conditions for decoupling of equations will be expressed in terms of these matrices. 



3 Decoupling of correlation functions 

The simplest correlation functions are one-point functions which determine the av- 
erage densities of particles of each type at each site. Thus we are interested in the 
equation of motion of the one-point function {rii(x)) := (8 T t x \i). This is in fact the 
probability that site x is occupied with a particle of type i. In the operator formalism 
this one-point function is written as a matrix element of the corresponding operator, 
namely: 

(m{x,t)) = (S\Ei{x)\P(t)) =: (Ei{x)), (21) 

where by the last equality we have defined the bracket or the average of an operator. 
From equation (||), we find 

Uni(x)) = m{x),H}) = {{E t {x)M^ ~ ^x)]) + {[Ei{x)Mx,x + I)}). (22) 
at 

Each of these two terms leads in general to two-point functions on their relevant sites, 
the first term on (x — l,x) and the second on (x,x + 1). We now ask under what 
condition a cancellation occurs in each of these terms separately so that we end up 
with only one point functions on the right hand site. A little reflection on the words 
written in italic shows that the question we are asking and henceforth its answer, 
does not depend at all on the underlying lattice and the range of the interaction, in 
so far as we are only considering Hamiltonians with two body interactions. 
To find the answer to our question we calculate one of these terms say the second 
one. We note that: 

{S\[E i {x),h(x,x + l)]=R%"{S\[E ii (x),E tia {x)E 1/f3 (x + l)} (23) 
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where we have ignored the diagonal part of h(x,x + l) which obviously commute with 
Ea(x). Using the commutation relations ( |i~4| ) and also the property ( |l5| ) we find: 

(Sl[i^(x),%,x+l)]=£i^ (24) 

The second term vanishes in view of equation ( |18| ) and we are left with: 

(S\[E u ,h(x,x + l)} = (R i ) a>p (S\E a (x)E p (x + l) (25) 

Expanding the right hand terms we find 

{S\[E u ,h(x,x + l)} = R\ {S\E (x)E (x + l)+B} j>0 {S\E j (x)E (x + l) 

+ R\ k {S\E (x)E k {x + 1) + K i j:k {S\E J (x)E k (x + 1), (26) 

where now we are using the summation convention on Latin indices and R-'^fc stands 
for (R 1 )^. We now use the operator identity 

E = 1 - Ex - E 2 E p (27) 

and eliminate the operators Eq(x) and Eq(x + 1) on the right hand side of the above 
equation and demand that all the quadratic terms vanish. It is easily seen that this 
cancellation occurs when the matrices R 1 satisfy the relations: 

^j,k = R-Vfc + ^\7>o ~~ R-'o.o V i, j, k. (28) 

This simply means that for each matrix R 1 , all the elements are fixed once the elements 
of the first row and column are determined. 



Similar calculations for the first term of (22) leads to the following condition: 

9 l j,k = SV + si i,o " SVo V i, j, k. (29) 



where the matrices S 1 have already been defined in (IS). 

Multiplying both sides of fl26|) by \P(t)) it is seen that once the above conditions 
are satisfied, the equations of motion of 1— point functions depend only on one 
point functions. A short calculation shows that once the equations of motion for one 
point functions are closed, it guarantees that the equation of higher order correlation 
functions are also closed, that is their equation of motion depend solely on the Appoint 
functions with k < n. In this way the hierarchy of equations of n— point functions is 
terminated and closed at any level and the system amends itself to exact analytical 
treatment. Every model whose rates satisfy the relations (28, 29) is solvable in the 
above sense. 

Moreover we show that the equations of motion of two-point functions at different 
times are also closed under the above conditions. In some models like voting models 
p7| ] the calculation of such correlation functions are important. 

To see this consider the two point two time correlation function (rii(x, t)rij(x' , 0)). We 
have: 

^(n i (x,t)n i (x , ,0)) = (S\e- tH [E u (x),H]e tH E n (x')\P(0)) (30) 
= (5|[^(x),^]e tH ^(x')|P(0)) (31) 
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Under the conditions (28 and |29|) , we know that the commutator in the above equation 
is expressible as the sum of local site operators. Restoring the operator e~ tH , i.e., 
(S\ — ► (S\e~ tH , we find that on the right hand side of this equation only two point 
functions appear. Therefore the equation of motion of two-time two-point functions 
will also be closed. 
Remarks: 

1: The emergence of two sets of constraint for the rates does have nothing to do with 
a given site x, having two neighbors x — 1 and x + 1 and hence on the underlying 
lattice. On a general lattice we still have only these two equations. It only reflects 
the fact that the interaction Hamiltonians defined on each link may not be symmetric 
under the interchange of its ends. Therefore each site x contributes to two types of 
interaction Hamiltonians, namely those of the type h(x,y) which for all different y's 
lead to ( |2"S| ) and those of the type h(z, x) which for all different z's lead to (p9|). In fact 
for symmetric models in which there is no driving force one expects that these two 
sets of constraint become identical. This is in fact the case, since for these symmetric 
models, one has: R^L = R V n" a and hence R 1 ^ = S 1 p ! a, which makes the two set of 
conditions identical. 

2: For general models the number of independent rates is (p+ 1) 4 — (p + l) 2 . To count 
the number of conditions on the rates we note that a matrix such as R 1 imposes 
p 2 linear equations on the rates, since its first row and column determine all the 
other elements. There are 2p such matrices and hence the number of conditions is 
2p 3 . Therefore the number of parameters of the solvable family we are considering is 
N p = (p + l) 4 — (p + l) 2 — 2p 3 . For the one species case we find Ni = 16 — 4 — 2 = 10. 
For the two and three species models we have respectively 56 and 186 free parameters 
respectively. 

3: It appears that the family of integrable multi-species models live on a manifold 
of huge dimension and one is at complete ease to choose many models of his own 
choice for adapting them to any physical situation. However this is an illusion and 
as we will see, in choosing physically interesting models from this manifold one is 
much more restricted than it appears at first sight. The reason is that we are usually 
interested in simple subclasses of these models, ones in which we can set many of the 
irrelevant or uninteresting parameters equal to zero either on physical grounds or to 
make our analysis simple and transparent. However it often happens that once we 
set some of these parameters equal to zero, the set of admissible parameters collapses 
drastically so that we are left with totally uninteresting models or models which 
are actually equivalent to one species models. Geometrically the constraints are a 
set of hyperplanes which pass through the origin. The admissible rates lie in the 
intersection of all these hyperplanes. However all the rates are also constrained to be 
positive. It may happen that this intersection of hyperplanes, although a manifold 
of high dimension, may intersect the positive sector of the space of parameters in a 
very low dimensional submanifold. As an example consider a model with say 6 rates 
r i) r 2> r 3; r 4> r 5 an d r§ and one constraint r\ + r<i + r% + + r§ — rg = 0. The set of 
admissible rates is obviously a five dimensional plane. However if for some reason we 
are interested in those models in which the rate r§ is vanishing, then we are left with 
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the constraint r\ + ri + + + r§ = which forces all these remaining parameters to 
vanish. In the multi-species case this difficulty shows up more severely since we are 
dealing with a large number of hyperplanes. Thus it is nontrivial to find physically 
interesting multi-species models in the class discussed above. 

We will conclude this section with the final form of the equation of motion of one 
point functions. 

Collecting the remaining linear terms in the equation ( p6| ) and its counterpart, we 
find (with summation convention understood for Latin indices): 

^(ni(x)) = R'o.ofaoOc) +n (x + l) - 1) 

+ B} jt0 ( nj (x)) + K l 0J (nj(x + 1)) 

+ (x-»x-l,R-»S) (32) 

This equation is specific to a one dimensional lattice. In a general lattice it should 
be modified appropriately, the modification is however straightforward (see section 
4). Before going to the consideration of a two-species model, it is instructive to 
recapitulate the findings for the one species model from this general perspective. 



3.1 The one species separable models 



If there is only one species of particles on the lattice, we have only two matrices 
namely R 1 and S 1 . These are two by two matrices, subject to the conditions: 



R 1 l,i + R 1 o,o — R 1 i,o + R 1 o,i 

S 1 !,! + S 1 o,o = S 1,0 + S 1 o,l 



(33) 
(34) 



When expanded by using the definition (28, |2^ ) of the matrices R 1 and S 1 , they 
yield: 



Ri\ + R[{ + Rqq + R t 



01 



00 

>11 

00 



d10 



+ -RqJ + R 



10 I oil 
10 ' - n io 



R 01 i nil i pOl , pll 
^01 + ^01 + ^10 + -"io 



(35) 
(36) 



One can now eliminate the diagonal terms from normalization to transform the above 
relations to: 



1.0 + ^10 + -^oo + -^00 



^00 , D 01 
pUU I pin I pin 

-"■01 - ft oi -"oo 



R 



10 



L R 11 _L "00 , p 01 

oi + Rq\ + K n + K n 



-"oo 



R 00 i plO , pOl 
-Kll + K ll + K 1Q 



R 10 



(37) 
(38) 



which are the relations given in [21]. In the next section we will consider a two species 
model. 



4 A two species separable model 

As mentioned above the rates of the solvable two species models of the type considered 
in this paper live on a 56 dimensional hyperplane. By exploring various regions of this 
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plane, one can find interesting two species models suitable for various applications. 
In this section we explore a small region of this plane by imposing extra symmetry 
requirements on the model. 

We consider systems in which there is no driving force, i.e. systems which have the 
symmetry R^ V a = R'g 1 ^- As noted above (see equations (28,^9), for these models 
the R matrices and the S matrices lead to identical constraints. Furthermore we 
restrict ourselves to those models in which an individual species is neither created 
nor annihilated but only changes its label. Such models may be appropriate for 
description of voting processes or the spreading of epidemics. This means that we are 
setting Rq° = R®'® = R^j = Rf Q = In other words the number of Latin indices 

2 

should be equal as subscripts and superscripts of R. A rate such as i? 1 ' means that a 
voter with vote 1 spontaneously (or due to the effect of environment) changes his or 
her vote to vote 2. (or a healthy individual 1 gets infected due to the interaction with 
the environment.) In a time interval dt two voters with different votes 1 and 2, pass 
each other, without changing their votes, with probability Rf^dt. It may also happen 
that on this close contact the voter 1 changes his or her idea and switch to vote 2. 
This will happen with probability R^dt. Note that the voting processes that have 
been studied in the literature contain only two species + and — with no vacant site. 
Here we have also vacant sites, and the voters can move in free space and interact 
with each other. 
Equations (|28|) now yields: 



plO 


+ -Kqi 


= R\\ 


+ R\i 


r>20 
K 10 


i R 20 
+ -^01 


— Z? 21 

- K n 


+ Ru 


(39) 


plO 


_i_ plO 


= R 12 


+ R\l 


r>20 

-fMo 


+ K 02 


- K 12 


+ Ril 


(40) 


n 2Q 


+ 


= R 21 


+ R% 


c>20 
n 20 


, R 20 




+ Rfi 


(41) 


plO 
n 2Q 


+ ^02 


— -"22 


+ R-22 


p20 
a 20 


i R 20 
+ K 02 


— P 21 

- -n-22 


+ Rll 


(42) 



The diagonal terms like R^, R^ ■ ■ ■ are a source of trouble, since they are minus the 
sum of a number of rates and we would better get rid of them. To do so we proceed 
as follows: We note that in each pair of equations above, the sum of the right hand 
side terms adds up to zero, due to normalization. Thus the sum of the left hand sides 
must also add up to zero. If we do so and substitute for the diagonal terms their value 
from the normalization (i.e: substitute R\q with —(R^q + + i?f{]), we find that 
the first and the last pair of equations of the above set lead to trivial identities. We 
also find that the second and the third pairs lead to one single identity, meaning that 
in each pair one is redundant. Thus we keep this last identity and safely ignore all 
the equations which contain diagonal terms. Therefore we are actually dealing with 
5 independent equations relating positive rates, the final forms of which are: 



p 20 , p 02 _ p 21 , p 22 plO i pOl _ pl2 , pll ( A o\ 

K 10 + K 10 - K U + K U K 20 + K 20 - K 22 + K 22 

p 20 , p20 _ p21 , p22 plO i plO _ R 12 , pll ( AA\ 

K 10 + K 02 - K 12 + K 12 K 20 + K 01 - K 21 + K 21 ( 44 J 

-^01 + -^01 = -^02 + -^02 ( 45 ) 

These are the final conditions on the rates for this kind of two species model. 
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4.1 Equations of motion 

To obtain the equations of motion for the above two species model, we use ( |32| ) and 
obtain for a d dimensional rectangular lattice with unit vectors e r ;r = 1, • • • d and 
with the abbreviation (rij(x)} — > rij(x): 



dt 
and 

dt 



d m(x) = R$V 2 n 1 ( X )-2d(R™ + R 2 l ° 1 )n 1 (x) 



+ ^V 2 n 2 (x)+2d(^+i?^)n 2 (x), (46) 



-|n 2 (x) = i^°V 2 n 2 (x)-2 c /(i?^ + i?^)n 2 (x) 

+ ^?V 2 n 1 (x) + 2d(i?[] 2 + J R^)n 1 (x). (47) 



where V 2 n(x) := (X^=i n ( x + e »*) + n ( x ~~ e ^) ~~ 2dn(x)j stands for the discrete <i 
dimensional Laplacian. 

All the terms in these equations can be understood intuitively. For example a term 
like i?Q 2 V 2 n 2 (x) measures the diffusion of particles of type 2 which change their type 
or color to 1 as they hop. 

If one begins to write the equations intuitively taking all the complex interactions 
into account, one finds many other terms. But at the end they all cancel out. For 
the description of the solution, it is convenient to define new parameters: 



71 
7 
D 

D' 



2d(i? 2 ? + <) 72 = 2^ + (48) 

7i + 72 (49) 

^01 + ^01 = ^02 + -^02 - 1 ( 50 ) 

dIO olO _ p20 7?20 /ki \ 

-"-01 - ft 02 — -"02 -"01) \ 0L ) 



where we have used the last relation of (43) in the last two relations and have rescaled 



time to set D = 1. These parameters have obvious physical interpretations, 71 and 
7 2 respectively determine the overall tendency of particles of type 1 and 2 to switch 
their types. D is the diffusion constant of the particles when we ignore their types and 
finally D' a kind of relative diffusion constant. It measures the difference of diffusion 
constant for the particles that do not change their type in hopping, compared with 
those which do so. 

Solving these coupled system of differential-difference equations will give the distri- 
bution of both types of particles in space and time. In principle it is possible to go 
to fourier space and diagonalize the resulting matrix equation. However it is better 
to proceed in a more physically transparent way by defining new densities: 

n(x) = ni(x) + n 2 (x) (52) 
0(x) = 7irai(x) - 7 2 n 2 (x) (53) 

Here n(x) is the total density of particles (when we ignore their types or colors) and 
0(x) is a weighted difference of densities. 
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It is now a matter of simple algebra to use equations (E^) to arrive at the following 
equations for these new densities: 

J^n(x) = V 2 n(x) (54) 

which means that if we ignore the color of particles, they perform simple diffusion 
with diffusion constant D. We also obtain 



J^(x) = iyv 2 0(x) + £>"V 2 n(x) - 70(x) (55) 

where 

D" := Rlhi ~ <72 - RX ~ « (56) 

The total number of particles of each species obey very simple equations which are 
obtained by summing the above equations over x. Denoting the number of particles 
of species 1 and 2 by N% and N 2 respectively, and noting that N 2 = N — N\ where 
N is the total number of particles, we find from (^5[) by summing (j) over all the sites: 



—Ni = -7iiVi + 72^2 = -71-^1 + 72(N - JVi) (57) 
the solution of which is: 

NAt) = 72 N + e~ (7l+72)f (iVi(0) — — N) (58) 

7i + 72 7i + 72 

N 2 (t) = 71 N + e~ (7l+72)f (iV2(0) — N) (59) 

7i + 72 7i + 72 

We now discuss the solution of equations (|54-55) which determine the spatial distri- 
bution of particles in time. 
Let us define the fourier transforms 

n(q, t) := ^ e iq - x n(x) 0(q, t) := ^ e iq ' x 0(x) (60) 

X X 

where q = (qi, q 2 , * * ■ 9d) an d for all i, % € [0, 2ir), from which we find: 

n(x) = n(q,t)^ *(x) = j^e-^^t)-^ (61) 

With the definition: 

r=d 

Q = 2j2(cosq r - 1) (62) 
r=l 

we find the equations of motion for these generating functions as: 

n = Qn (63) 
= Q(D'0 + D"fi) - 70 (64) 
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with the general solution: 



n(q,t) = n(q,0)e Qt (65) 
ftq,t) = e «^)*[^( g ,0)-_^_fi(g > 0)] 

Once the initial distributions of particles of each type is known, these equations al- 
low us to determine the distributions of both types of particles in later times. If the 
particles are in a finite volume, then the above solutions are still valid, except that 
the momenta q will take discrete values. 

The large scale behaviour of these densities is determined by going to the limit of 
q — > where Q — ► — |q| 2 . For illustration we consider two simple examples. 

Example 1: Let us assume that the particles (voters) change their type (votes) 
only on encounter with other particles (voters). In this case we have: R^q = it^o = 
^20 = ^20 = °- From (H) we find that the onl y 

non-zero parameters which remain 
are R% =: P, R\l = R\\ =: A and R% = R% =: D, subject to a relation: D = P + A, 
where we have introduced new simple labels for the rates. The labels D, P and A 
stand respectively for " Diffusion ", " Pass" and " Agreement": 

(67) 
(68) 
(69) 
(70) 
(71) 

In this case we find from (46 and 47) that both types of particles diffuse through each 



1 





<-» 


1 


with 


rate 


D 


2 





«-► 


2 


with 


rate 


D 


1 


2 


2 


1 


with 


rate 


P 


1 


2 


-» 2 


2 


with 


rate 


A 


1 


2 


-» 1 


1 


with 


rate 


A 



other without any interaction: 

-ni(x) = D\7 2 ni (x), dt n2 ^ = DV2n2 ^ ^ 

However this is peculiar to one point functions and the equations of two point func- 
tions will indeed be coupled to each other by interaction parameters. This is an 



example of an observation first made in [32| according to which some hamiltonians 



may lead to the same set of equations for one point functions. Here our hamiltonian 
is equivalent to a free hamiltonian as far as the one point functions are concerned. 



Example 2: Let us now assume that one of the particles say type 1 does not change 
its type, or one kind of voters is persistent in his vote, that is, = ii^o = 0- Fne 
first equation of ( |43"|) then leads to Rl\ = = 0. In this case we find from d48| ) that 
D" = and D' = 1. The relations (^3|) reduce to the following relations between the 
remaining rates of reactions: 

p20 r>21 i r>22 /''ZQ'v 
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^01 = ^02 + ^02 ( 74 ) 
-^20 + "^01 = R 21 + ^21 ( 75 ) 
-^20 + "^20 = -^22 + -^22 ( 76 ) 

For simplicity we first consider the large scale form of the distribution functions. From 
(|65|) we obtain 

n(q,t) = n(q,0)e-l q l 2 ' (77) 
0(q,i) = 0(q,O)e(-l q l 2 -^ (78) 

Let us consider a situation at time t = where there are N 2 particles of type 2 at the 
origin, i.e. ?i2(x, 0) = N 2 5(x) in a uniform see of particles of type 1 of density p, i.e. 
?ii(x, 0) = p. Thus the initial values of fourier transforms are: 

ni(q, 0) = p(2vr) d «5(q) n 2 (q, 0) = N 2 (79) 



or 



n(q, 0) = N 2 + p(2Tr) d 5{ci) 0(q, 0) = - 72 iV 2 . (80) 
Inserting these into (0), we find: 

n(q,t) = (N 2 + p(2ir) d 5(q))e-^ 2t = N 2 e-^ 2t + p(2ir) d 5(q) (81) 

i) = - 72 AT 2e -(72+|qlX ( 82 ) 



Taking the inverse fourier transform and using the definition ( ]52| ) we find: 

ni(x,t) = p+-^- r e^(l-e-^) (83) 
(4vrt)2 

iV 9 -|x| 2 

n 2 (x,i) = -^e— «r»« (84) 
(4?rt)2 

(85) 

The initial number of particles of type 2 diffuse and gradually the particles of type 2 
turn into type 1. At the end, no particle of type 2 remains. 

We can also consider the small scale behavior of these distributions. In this case the 
solution of the equations ( |63| ) is given by modified Bessel functions. We give below 
the solution for the initial distribution of one particle of type 2 at the origin and a 
uniform distribution of particles of type 1 at all other sites except the origin. That 
is ( n 2 (x,0) = <5 Xj0 or n 2 (q,0) = 1) and ( m(x) = p(l - £ Xj0 ) _or ,(ni(q,0) = 
p((2ir) d 5(q) — 1). For these initial conditions we can find n(q, t) and c/>(q, i) from (|65]). 
The result is: 

n(q,t) = (1 - p)S l + (2tt) V(q) (86) 
0(q,t) = -^e-W-w)* (87) 
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On the lattice we will find: 



n(x,t) = p + (l-p)^^V-^ q 

These integrals can be written in terms of Modified Bessel functions 
I x (t) = ± e-^ x+cosqt dq. Thus we find: 

= /) + (l-p-e-^)e^n^x(2t) (90) 

X 

n 2 (x,i) = e-^e- M *n^x(2t) (91) 



5 Discussion 

We have obtained the general condition under which the hierarchy of equations for 
n— point functions of a multi-species reaction diffusion system truncates at all orders, 
that is the equations of motion of n— point functions depend only on those of lower 
correlation functions. Also under the above conditions the equations of motion of 
two-point functions at different times form a closed system. We have selected out a 
class of 2 species models which may be appropriate for the description of the spread 
of an epidemic or as a voter model. This work can be extended in several directions, 
even if one restricts oneself to the two species models. First, one can study other 
two species models which include coagulation, de-coagulation or birth and death 
processes. Second, in the two species model one can change our interpretation of 
particle as a hole and take it to be a real particle. In this way one can adapt the rates 
for description of other models. Third, one can study also the asymmetric models 
and relax our simplifying assumption on the symmetry of rates, and finally one can 



consider the similarity transformation [28, 29, 3C, 31] on the model to obtain other 
exactly solvable models. 
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